/*!  \file  FormSVD.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      09/29/2013 01:39:40 PM
 *   \copyright GNU Public License.
 */
#include <iostream>
#include <assert.h>
#include <stdlib.h>
#include <iomanip>
#include <string.h> // memcpy
#include <math.h>
#include <stdio.h>

#include "Matrix.h"
#include "FormSVD.h"
#include "BiSVD.h"
#include "BiReduction.h"
#include "commonSVD.h"
#include "setThreads.h"
using namespace std;

FormSVD::FormSVD(const char * name_) : Implementation(name_, FormSVDGroup) {
	setCompileTime(__DATE__, __TIME__);
	canBeParallel = true;
};
bool FormSVD::checkInput(Matrix & Q, Matrix & P, Matrix & U, Matrix & V) const{
	bool res = true;
	res &= Q.isReflector() && Q.isSquare();
	res &= P.isReflector() && P.isSquare();

	res &= V.isOrthogonal() && V.isSquare();
	res &= U.isOrthogonal() && U.isSquare();
	
	res &= Q.rows() == U.rows();
	res &= P.rows() == V.rows();
	// A = Q B P^T
	// B = U S V^T
	// A = (Q U) S (P V)^T
	return res;
}
bool FormSVD::setNumberOfThreads(const int nThreads_) {
	assert(canBeParallel);
	SET_NUM_THREADS(nThreads_);
	nThreads = nThreads;
	return true;
}
// OR = orthogonal
// G  =  generate
// BR = from Bidiagonal reduction 
ORGBR::ORGBR(): FormSVD("Lapack orgbr + dgemm") {
	setCompileTime(__DATE__, __TIME__);
}
ORMBR::ORMBR(): FormSVD("Lapack ormbr") {
	setCompileTime(__DATE__, __TIME__);
}
void ORGBR::compute( Matrix & Q, Matrix & P, Matrix & U, Matrix & V ){
	if (!checkInput(Q,P,U,V) ) return;
	Q.generateOrthogonal();
	P.generateOrthogonal();
	U = Q * U;
	V = V * P.T() ;
}
void ORMBR::compute( Matrix & Q, Matrix & P, Matrix & U, Matrix & V ) {
	U.applyReflectorFromLeft(Q);
	V.applyReflectorFromRight(P.T());
}
FormSVDUnitTest::FormSVDUnitTest(FormSVD & form_): form(form_) {
	int l = strlen(form.getName());
	char buf[l + 20]; 
	sprintf(buf, "FormSVD %s", form.getName());
	setName 				(buf);
	setDescription  	("Test registered bi form implementations");
	setTestAuthor		("Yordan Madzhunkov");
	setDeveloperAuthor("Yordan Madzhunkov");
	setCompileTime(__DATE__, __TIME__); // each class that inherits this one should also have this line.
	setCompileTime(form.getCompileTime());
};
bool FormSVDUnitTest::test() const { 
	ImplementationFactory * factory = ImplementationFactory::instance();

	Implementation * imp = factory->first(BiReductionGroup);
	if (!imp) return false;
	BiReduction & reduction = dynamic_cast<BiReduction&>(*imp);

	imp = factory->first(BiSVDGroup);
	if (!imp) return false;
	BiSVD & bisvd           = dynamic_cast<BiSVD&>(*imp);

	const int m = 5, n = 6;
	double a[] = {8.79,  6.11, -9.15,  9.57, -3.49,  9.84,
	              9.93,  6.91, -7.93,  1.64,  4.02,  0.15,
	              9.83,  5.04,  4.86,  8.83,  9.80, -8.99,
	              5.45, -0.27,  4.85,  0.74, 10.00, -6.02,
	              3.16,  7.98,  3.01,  5.80,  4.27, -5.31};

	Matrix A   (n, m, General, a); 
	Matrix B   (A.rows(), A.cols(), BiDiagonal );
	Matrix Q   (A.rows(), A.rows(), ElementaryReflectors);
	Matrix P   (A.cols(), A.cols(), ElementaryReflectors);
	Matrix tmp (A, CopyWithData);
	reduction.compute(tmp, Q, B, P ); // tmp is destroyed in this computation
   Matrix S(B.rows(), B.cols(), Diagonal);
	Matrix U(B.rows(), B.rows(), Orthogonal);
	Matrix V(B.cols(), B.cols(), Orthogonal);
	bisvd.saveCompute(B, U, S, V);
	form.compute(Q, P, U, V);

	const int  nNorms = 6;
	double norm[nNorms];
	const double limit = 1e-12;
	getNorms(norm, A, U, S, V);
	bool res = true; 
	for(int i = 0; i < nNorms; i++) {
		res &= norm[i] < limit;
	}
	return res;
}
REGISTER(ORMBR);
//REGISTER(ORGBR);
